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O ; Abstract 

In this paper I describe a specialized algorithm for anisotropic diffusion determined by a field 

Q-f 

of transition rates. The algorithm can be used to describe some interesting forms of diffusion that 
occur in the study of proton motion in a network of hydrogen bonds. The algorithm produces data 
that require a nonstandard method of spectral analysis which is also developed here. Finally, I 

O" 

^ ■ apply the algorithm to a simple specific example, 
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I. INTRODUCTION 

Protons migrating in water have an anomalously high mobility [l| and their diffusion is 
actually limited by the continuous rearrangement of hydrogen bonds . Indeed protons 
migrating in ice move faster than protons in water, as the transition rate from one water 
molecule to the next is enhanced by the higher molecular order in ice. Proton mobility 
increases whenever water molecules are constrained, as in carbon nanotubes jj, 4]. Local 
electric fields also orientate water molecules, and thus should lead to a local increase of proton 
mobility, and indeed it is now known that there is a definite water dipole orientational order 
in the hydration water close to ionizable residues in hydrated proteins : this is 

a collective property, which is somewhat independent of the individual fluctuations of the 
water dipoles. 

Here I am not concerned with the detailed simulation of proton motion which is the sub- 



ject of several specialized papers like 
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but rather I wish to set up the framework 



for a simulation of the random walks performed by protons in some interesting context, like 



proton migration on the surface of hydrated proteins [13|. The basic idea is that protons 
move faster in the network of hydrogen bonds just where there is a higher molecular order, 
i.e., the transition rate is higher where there is higher spatial order, and because of the 
continuous rearrangement of the water molecules which make up a fluctuating bond struc- 
ture, the random walk performed by a single proton can actually be viewed as a walk in 
continuous space and continuous time, as long as the time resolution of the process is longer 
than the relaxation time of water dipole motion. 

Here I take for granted that there is some induced order in the hydrogen bond network, 
like the dipole field described in , and I introduce a corresponding field of transition 
rates 7(1*, i), such that the time-dependent probability density p(r,t) and the associated 
probability Ap = pAV of finding a random walker (a proton) in the small volume AV at 
position r and time t, yield the following equation for the decrease of Ap, due to random 
walker escape from the region, 



-7(r,t)Ap(r,«) (1) 



I also assume that 7(1", t) is a continuous, differentiate function. 

The situation is illustrated in figure [TJ which shows a random subdivision of a plane 



2 



region: a set of positions - marked by the large black dots - is associated to small surrounding 
regions; the arrows in the figure mark the flow of random walkers in the central region to and 
from the bordering regions. If the area (actually, the line length in this 2D representation) of 
the interface between the central region and the kj-th region is A^. , and the total interface 
area of the kj-th region is A\. , then it is easy to see that the total derivative of Api is 

dApi v^A i)fc 

3 3 

and the global flow in this discretized system is described by a system of coupled linear 
differential equations. 

More importantly, we can define currents for the inflow and outflow of random walkers 
from a modified form of Fick's law 

j liPi ~ IkPk fo s 

Ji^ k = a (3) 

where A^k denotes the distance between the centroids of the bordering i-th and k-th region, 
and the parameter a is akin to the diffusion coefficient, but is measured in different units 



(it has the dimensions of a surface) [l4]. The current in the previous formula is actually a 
projection along the direction that connects the centroids of the bordering region and it is 
easy to generalize to the continuum case and find the outflowing current 

J = -«V( 7 p) (4) 

so that finally one finds the following Fokker-Planck equation from the conservation of the 
total number of random walkers 

% = «V 2 ( 7 p) (5) 

assuming that a does not depend on position. 

In the following sections I describe an algorithm to simulate this kind of diffusive motion: 
first I discuss the angular distribution, then confinement to motion on surfaces, and in 
section IIVI I show how to extend the algorithm for asynchronous updates. In section [V] I 
give a recipe to analyze asynchronous data. In section IVII I discuss a simple example, and 
finally in section \VU\ I give a short summary and outlook for the utilization of the algorithm. 
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II. ANGULAR DISTRIBUTION 



From equation (T4J) we see that the current actually contains two contributions 



J 



aV(7p) 



a (pV7 + 7V p)) 



(6) 



however when we consider the problem at hand - namely, the diffusion of protons in the 
network of hydrogen bonds, and we remark that we wish to describe the individual proton 
motion, then we notice that we are only interested in situations where Vp = 0. In fact, pro- 
tons repel other protons that are too close, and obey a sort of effective exclusion principle - 
which is actually independent of their fermionic nature - and the position of the individual 
proton corresponds to a peak of the instantaneous probability density: therefore the current 
defined in (J3J) has the same direction as V7 in all cases of practical interest. This direction 
corresponds to the average proton motion, but for a single transition to a nearby site it can 
only define the axis of an angular probability distribution. Here I make the simplest possi- 
ble choice, namely that the angular probability distribution is a simple dipole distribution 
defined by the normalized conditional probability density for the unit vector n 



where Jo is the normalization factor (Jo = 2tt in the 2D case and Jo = 4-7T in the 3D case), 
and n is the unit vector 



so that the decrease of the density p due to the flow in the angular range AQ, during a given 
time interval At, is 



The constant inside the parenthesis corresponds to the isotropic loss term, while the other 
term is associated to the current (J6]). From a comparison of the elementary flows of random 
walkers in direction n we find 




(7) 




(8) 




(9) 



J • n = 7p 



AP 



(10) 



n • n 



so that 



AP = I a 



|V 7 | 



(11) 
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and the conditional angular probability density is 

P(n|n ) = J o ( K 1 + J o«^ n ■ n o) (12) 

The conditional angular probability density (fT2|) can be used to generate random walks 
discarding the time information. Here I take the following time-independent expression for 
the transition rate 

7M)=7(r) = p^ + 5 (13) 

which has an obvious symmetry center, located in the origin, which corresponds to the po- 
sition of the peak value as well. This transition rate is motivated by the considerations put 
forward in the introduction: if protons migrate in a hydrogen bonded network with polariza- 
tion centers that create partial ice-like order in their neighborhood, then the transition rate 
(TTBl is highest, and saturates, close to the polarization centers, and decays to a constant 
value with a 1/r 2 behavior for r ^> T (i.e., it has a radial dependence like the potential of 
electric dipole fields). Notice also that the anisotropy coefficient is 

IVtI = (14) 
7 |if + r 2 1 ; 

The techniques to generate random angles which are distributed according to the probability 
density (|T2|) are reviewed in appendix and figures [2]l5] show some examples: in these 
examples all length and distance units are in arbitrary units. Figure [2] shows random walks 
around a single center with transition rate (fT3|) : the random walker starts at the origin, 
with a fixed step length = 0.005 arbitrary units; the horizontal and vertical scales are also 
labeled with the same arbitrary length units; the parameters of the transition rate function 
are the same in these simulations, A = 1, B = 0.1, and T = 1, while a changes in the three 
cases displayed in the figure. Larger values of a correspond to higher anisotropy, and we see 
that as the anisotropy grows, the random walk becomes more and more compact. 

Figure [3] shows a random walk with two centers at positions ri = (—1,0), r 2 = (1,0) 
(arbitrary units): the random walker starts at the origin, with a fixed step length = 0.01 
arbitrary units; the horizontal and vertical scales are also labeled with the same arbitrary 
length units. In this case the transition rate is similar to ( fT3l) . but with two centers, 



A 1 A 



7(rH |r _ ri| ' 2 + r; + |r ^ + ri +B (15) 
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with A\ = A% = 1, B = 0.1, and I\ = r 2 = 1, and a = 1. Here the random walker explores 
the regions around both centers. 

Figure H] shows a situation which is similar to figure [3J although it is more complex. The 
transition rate is once again similar to (IT5|) . but now it has ten centers, 



with Ak = 1, B = 0.1, and = 0.1, and a = 0.025; the step length is = 0.01 arbitrary 
units. The centers are scattered randomly, with a lower bound on the minimum distance 
between them; the figure shows three snapshots at different times in the simulation, as the 
random walker starts from the center of the figure, drifts to one of the centers and later 
migrates to other neighboring centers. 

Finally figure [5] shows a random walk in space about two centers at r x = (—1,0,0), 
Y2 = (1,0,0) (arbitrary units), which is very similar to the random walk in figure [3j the 
transition rate is still given by expression (fT5l) . with A\ — A 2 — 1, B — 0.1, and Ti = T 2 = 1, 
and a = 0.5, with a fixed step length = 0.1 arbitrary units. Once again the random walker 
explores the regions around both centers. 

III. DIFFUSION ON SURFACES 

In many cases it is important to confine the motion of the random walkers to some 
particular portion of space, for instance in the case of protons on hydrated proteins the 
motion is confined to the thin hydration layer. The simulation method outlined in the 
previous section can be adapted to provide such a confinement to a surface: in this case one 
can define at each step the tangent plane at the position of the random walker and proceed 
as in the 2D case. Obviously the gradient of the transition rate 7 used in the formulas of the 
previous section must be projected on the tangent plane, and moreover the directions must 
be generated according to the 2D angular distribution (see appendix [B]) . Figure [6] shows 
a random walk on a spherical surface with 10 centers as in the examples in the previous 
section: the sphere has radius 1 (arb. units), the transition rate is given by expression ffT6l) . 
with Ak = 10, B = 1, and = 0.25, and a = 0.2; the step length is = 0.1 arbitrary units. 




(16) 



6 



IV. ASYNCHRONOUS TRANSITIONS 

In the previous sections we have discussed the space behavior of the random walks, but 
obviously we can use the transition rate function 7(1*, t) to describe the time behavior as 
well. It is possible to choose a fixed time step At and use the transition rates to compute the 
probability of generating a transition in the time interval (t,t + At) (synchronous update), 
however this is inconvenient if the function 7(1", £) spans a wide range of values, because it 
means that the choice At -C min(l/7(r, t)) which is required for an accurate simulation, 
produces very long waiting times where the transition rate is very small (and therefore, very 
large amounts of sampled data). It is actually much more practical to use the transition 
rate to generate directly the transition times of each step, which we assume to be indepen- 
dent (asynchronous update) from one another. With this - rather natural - assumption of 
independency, it is very simple to generate the transition times, as explained in appendix iBl 
although this leads to uneven sampling, and requires a specialized form of spectral analysis. 

V. FOURIER ANALYSIS OF ASYNCHRONOUS DATA 

Using asynchronous sampling times it is not possible to use the standard Fourier or other 
similar spectral analysis techniques [ijj]. However the signal produced by the time-domain 
simulation is "exact", at least in the sense that there are no algorithmic artifacts due to 
sampling and it is desirable to extract as much information as possible. To this end, I 
notice that any function of the position of the random walkers must be stationary between 
successive transitions, and that it is possible to make direct use of the definition of Fourier 
transform 



where s(t) is any signal produced in the simulation, which depends on the positions of the 
random walkers (e.g., a component of the electric dipole moment if the random walkers are 
charged particles). The signal s(t) has the fixed value s n in the time interval (t n ,t n + 1), 




(17) 
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where t n is the time of the n-th transition, and we find: 




■tn+l 



n 



sin [u(t n+1 - t n )/2] 

U\tn+1 — t n )/ 2 



(18) 
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(19) 
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Using equation (fl9j) the Fourier transform can be evaluated exactly for all frequencies, and 
without aliasing: in practice this is possible, practical, and actually useful only for a small 
finite set of frequencies. If we had used a Discrete Fourier Transform (DFT) algorithm 
l^l l to analyze N real samples, we would have found N/2 independent Fourier coefficients, 
and using a Fast Fourier Transform (FFT) algorithm the time complexity of the calculation 
would be O(NlogN). If we use the algorithm defined by equation ([19]) to compute M values 
(M < N/2) of the Fourier transform, the time complexity is clearly O(NM), so that in a 
practical calculation we can only compute a reduced number of Fourier coefficients. However, 
I remark that in addition to being exact, the algorithm has another major advantage over 
the standard DFT calculations: there is no limitation to the set of frequencies that can be 
computed, and in particular one can choose a set of frequencies that is not evenly spaced and 
that is denser close to the origin, which is particularly useful in this case since the random 
walk - when considered as a noise process - is expected to produce a spectrum with a large 
power-law peak at low frequencies. 

This kind of analysis is actually limited by the finite time span of the generated signal: 
we see from equation ( Tl8|) that the Fourier transform of the generated signal is a sum 
of sine functions, and therefore it is not useful to represent the transform for frequencies 
lower than uo min = 71 /T where T is the signal duration and w m j„ is the lowest positive 
zero of the corresponding sine function. With this limitation, we can sample the Fourier 
transform at frequency values that are evenly spaced on a logarithmic scale and obtain a 
better representation of the transform close to the origin than is possible with conventional 
methods. 

In this approach we evaluate the Fourier transform of the simulated signal in the time 
interval (0, t^) and we implicitly assume that the signal vanishes outside this interval: this 
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is different from the standard (implicit) assumption in standard DFT analysis, where the 



observed signal is assumed to repeat periodically outside the observation interval |15j . If we 
introduce a rectangular window with a width equal to the observation interval (tjv)j we see 
that the present method returns a Fourier transform that is the convolution of the transform 
of the signal with the transform of the rectangular window, which is 

exp(-iut)dt = t N exp(-itut N /2)—± — (20) 

(ut at/2) 

As a consequence of the convolution associated to the rectangular window we see that a 
constant nonzero level produces a sharp peak centered at zero frequency, with a shape given 
by equation (1201) ; this peak corresponds to the standard DC peak in DFT analysis, and has 
tails with al/u 2 spectral behavior that may mimic the low- frequency behavior of a standard 
Debye relaxation with a very small decay rate. The mean level of the simulated signal is 

s = — y~) s n / dt = — y~] s n (t n+1 - t n ) (21) 

*N — Jt tN — 

and thus we can correct for the DC peak by subtracting its transform 

2s 

— exp(— lutx/ty sin(u;tAr/2) (22) 

UJ 

from the signal transform. 



VI. RANDOM WALK ABOUT A SINGLE CENTER 

As an example, I consider here a complete simulation (3D space and time data) for 
random walks about a single center (as defined by the transition rate ([TBI ) located in the 
origin. The transition rate function used in this example is given again by expression ( [TBI , 
with A = 100, B = 10- 8 , r = 0.1, a = 0.001, and with a step length = 0.01. The values 
of A and B have been chosen to maximize the range of 7 sampled by the random walker, 
while still keeping a rather short simulation time. I have generated 500 random walks and 
10000 transitions for each walk; the random walker always starts at the origin (where the 
center of ( fl3l) is located). The results of the simulation are shown in figures [TlfTTl figure [7] 
is the superposition of a few walks, and it is qualitatively clear that the density profile is 
very similar to that of the standard random walk in the plane. Figure [8] shows instead the x 
projection of the position signal vs. time, and this is not very different, e.g., from an electric 
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dipole component if the random walkers are charged particles. The insets show parts of the 
signal with increasing magnification, and the last inset displays clearly the stationary parts of 
the signal between successive transitions. Figure M shows the (unnormalized) distribution of 
time intervals between successive transitions: the figure demonstrates clearly that although 
the times have been generated according to the interval distribution in appendix [Bj the 
distribution in figure [9] is not a simple exponential, but rather it contains two different 
power-law regions (marked by the dotted lines in this log-log plot), which reflects the way in 
which the transition probability function is sampled by the random walks. One well-known 
property of ordinary random walks is that their mean square radius is proportional to time, 
i.e., (r 2 ) oc t: here we see (figure ITU]) that this linearity is recovered only asymptotically, as 
random walkers explore regions that are far away from the origin. Finally, I have used the 
power spectral estimation method of section |V] to analyze the x position signals (as those 
in figure [8]): the result is shown in figure [TT1 Figure [TTa is the power spectrum obtained in 
a single realization of the random walk, while figure [Tib is the average of 400 spectra. In 
each part of figure [UJ, the thin gray line represents an ideal power-law spectral density with 
the same slope as the the average of 400 spectra, i.e., a l// 2 spectrum, which is usually 
expected in these types of processes. In fact the random walks simulated here effectively 
sample asymptotically only a rather limited range of transition rates - even though the A/B 
ratio is quite high - and this means that the usual superposition argument that leads to 
more general power-law spectra does not apply here and it is quite natural to find a 
1/ f 2 spectrum. 



VII. DISCUSSION 



Before concluding this paper it is important to note that the correct continuum formu- 



lation of the diffusion equation in an inhomogeneous environment has been the subject of 
much discussion in the past and is still debated (see [l?], [3], and references therein). The 
generalization is unclear because the microscopic details seem to matter 171 ] . Moreover, 
there is also some interest towards the diffusion equation in various forms of anomalous 
diffusion 19j. Here I wish to stress again that the results presented in this paper are spe- 
cialized and are meant to address diffusion in structures like those described in 5, f|, unlike 
other approaches described in the existing literature that deal with more general diffusion 
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problems 120 



21 



221 ]; still, the diffusion equation (jSJ) is similar to equation (5) in reference 



171 ] . and therefore it is interesting to give one further look at its structure. In section HT1 we 
have seen that the current (pEJ) has two components, and the gradient term that generates 



the random walks discussed here rough 
(using the terminology of reference 



y corresponds to the so-called "spurious" drift term 



171]). The other term in the current, namely — ceyVp 



has been neglected because the single random walkers considered here are charged fermions 



pro- 



and obey a sort of effective exclusion principle. Indeed, in a context like that of 
tons repel because of their charge, while their spin structure, and therefore also their true 
fermionic character - and the Pauli exclusion principle - do not matter much; this effective 
exclusion principle has been used in the past for an Ising-like modeling of proton motion, 
where the presence or absence of a proton at a given position is treated like a pseudo-spin 
variable (see the discussion in the review paper 12j). The situation would be very different 
if space could be filled by a cloud of random walkers: in a case such as this - which roughly 
corresponds to random walkers that effectively behave as bosons - the neglected term should 
be included. Thus the actual importance of the different terms of the diffusion equation (]5]) 
depends on the bosonic or fermionic character of the random walkers. 

One prominently missing term in the diffusion equation (jSJ) is the usual drift term asso- 
ciated to external fields, however if we look at the structure of the current ([51) , we see that 
we can easily produce a flow unbalance with a space- (and possibly time-) dependent a, so 
that we obtain a modified diffusion equation 

dp 



Of 



V • [aV(7p)] = Va • V(7p) + aV%p) 



(23) 



with an additional a-dependent drift term. 

A final, important comment, is that the algorithm presented here has a sort of backward 
approach with respect to other existing algorithms for random walks and diffusion in inho- 
mogeneous environments, as it starts directly from the transition rates, instead of deriving 
them from a given diffusion equation (as, e.g., 20]). 

To conclude, in this paper I have described a novel algorithm for anisotropic diffusion, 
which is continuous both in space and time, and I have discussed its application to a simple 
example, in anticipation of further work that shall be carried out in a realistic simulation of 
noise in proton conduction in weakly hydrated proteins 23]. 
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APPENDIX A: ANGULAR DISTRIBUTIONS 



Here I consider first the planar angular distribution defined by the normalized probability 
density 

P(tp;g) = ^-(1 + gcostp) (Al) 

Using the standard inversion method (described in many textbooks, see, e.g., [2^] ) , one finds 
that the solution ip of the nonlinear equation 



1 / N 1 

— {<p + g sin <p) + - 



y 



(A2) 



has the distribution described by fl All) if y is a uniform variate on the (0, 1) interval. 
The generation of a random direction in space from the normalized probability density 

1 



P(n;g) = +#n- n ) 



(A3) 



requires two angles, a zenithal angle 9 and azimuthal angle ip, where n defines the zenithal 
axis, so that the probability of finding a unit vector n in the interval (9, 9+d9), and (<p, ip+dtp) 
is 



dP(9, ( p;g) 



1 



1 + g cos 9) d cos 9 



2ty 



dip 



(1 + gx) dx 



2ty 



dip 



(A4) 



i.e., the probability density is the product of two independent densities, one uniform with 
respect to <p on the (0,2-k) interval, and the other linear with respect to x — cos 9 on the 
(—1,1) interval. Using again the inversion method, one finds that 



2 

x — — 
9 




(A5) 



has the required linear distribution if y is a uniform variate on the (0, 1) interval. Since n is 
not usually parallel to the z direction, two rotations are also required: one first rotates the 
reference frame so that n is parallel to the z axis, this is followed by the angle generation 
step, and finally one must transform back to the original reference frame. 

APPENDIX B: TIME DISTRIBUTION 



Time transitions are generated according to the exponential distribution, which has the 
probability density 



dpt 
dt 



7exp(— jt) 



(Bl) 
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where 7 is the transition rate. The standard inversion method [24| can be used again to 
generate times 

t = —1ny (B2) 

7 

that are exponentially distributed if y is a uniform variate on the (0, 1) interval. 
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FIG. 1: Here space has been subdivided in elementary regions (the dots represent the centroids of 
the elementary regions): each region is characterized by a particular transition rate 7$ and by an 
elementary volume V*. In this figure the transitions between the central region i and the adjacent 
regions (denoted here by an index kj) are marked by the gray arrows. 
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FIG. 2: Random walks around a single center as explained in section [TIJ The random walker starts 
at the origin, with a fixed step length = 0.005 arbitrary units; the horizontal and vertical scales are 
also labeled with the same arbitrary length units. The width parameter is T = 1 in all simulations, 
and the other parameters have fixed values as well, A = 1, B = 0.1, while a. a = 1; b. a = VlO; c. 
a = 10. Larger values of a correspond to higher anisotropy, and we see here that as the anisotropy 
grows, the random walk becomes more and more compact. 
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FIG. 3: Random walk with two centers, as explained in the text. The walk is superposed on the 
contour plot of the transition rate 7(1"). The random walker starts at the origin, with a fixed step 
length = 0.01 arbitrary units; the horizontal and vertical scales are also labeled with the same 
arbitrary length units. 
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FIG. 4: Random walk with several (10) centers, a. after 10000 steps; b. after 40000 steps; c. after 
70000 steps. The walk is superposed on the contour plot of the transition rate 7(r). 
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FIG. 5: Random walk in space with two centers, located at the positions marked by the gray 
arrows. 




FIG. 6: Random walk on a spherical surface, with several (10) centers, some of which have not yet 
been visited by the random walker. 
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FIG. 7: Superposition of several random walks about a single center in the origin. Isodensity 
contour lines are superposed on the density plot. 
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FIG. 8: This figure shows the x projection of the position signal (this is not very different from 
an electric dipole component if the random walkers are charged particles) vs. time for one of the 
random walks in the example of section IVI( both position and time are in arbitrary units. The 
insets show parts of the signal with increasing magnification, and the last inset displays clearly the 
stationary parts of the signal between successive transitions. 
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FIG. 9: Unnormalized distribution of the time intervals At between transitions in the set of 400 
random walks of 10000 steps each described in the text: a. logarithm of the relative frequency 
vs. At (both in arbitrary units), which shows that the distribution is not a simple exponential, 
but rather contains two different power-law regions (dotted lines); b. log-log plot of the same 
distribution, where the two power laws are identified by the nearly straight sections. 
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FIG. 10: Mean squared distance (r 2 ) vs. time. In an ordinary random walk the mean squared 
distance is a linear function of time: here we see that linearity is recovered only asymptotically, as 
the random walkers explore regions that are further away from the origin, where the 7 is nearly 
constant. 
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FIG. 11: Spectral density calculated with the method of section [Vj sampled at logarithmically 
spaced frequencies: a. spectrum obtained from the x position signal for a single random walk; 
b. average of 400 spectra. Spectral densities and frequencies are in arbitrary units. The thin 
gray line represents an ideal l// 2 spectrum, which is expected for this kind of processes: the 
computed spectrum deviates from the ideal spectrum only at very low frequency, because of the 
limited observation time. Notice also that there is no upward bend at very high frequency - a hint 
of the absence of aliasing. In part a. it is clearly visible that the spectrum is sampled at (500) 
logarithmically spaced frequencies, because there is no crowding at the high end of the spectrum, 
and no rarefaction at low frequency, unlike spectra obtained with the Fast Fourier Transform or 
other similar algorithms. 
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